Ground state and dynamics of the biased dissipative two-state system: Beyond 

variational polaron theory 
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We propose a ground-state ansatz for the Ohmic spin-boson model that improves upon the varia- 
tional treatment of Silbey and Harris for biased systems in the scaling limit. In particular, it correctly 
captures the smooth crossover behaviour expected for the ground-state magnetisation when moving 
between the delocalised and localised regimes of the model, a feature that the variational treatment 
is unable to properly reproduce, while it also provides a lower ground-state energy estimate in the 
crossover region. We further demonstrate the validity of our intuitive ground-state by showing that 
it leads to predictions in excellent agreement with those derived from a non-perturbative Bethe- 
ansatz technique. Finally, recasting our ansatz in the form of a generalised polaron transformation, 
we are able to explore the dissipative two-state dynamics beyond weak system-environment coupling 
within an efficient time-local master equation formalism. 



I. INTRODUCTION 

The thermodynamics and non-equilibrium dynamics of 
quantum systems in contact with environmental degrees 
of freedom is a topic of primary importance in physics 
and chemistry^— and is also becoming increasingly rel- 
evant in biology for systems in which quantum effects 
may play an important role£~— In typical experiments, 
it is impossible to observe the degrees of freedom of the 
environment. The unmeasured correlations which build 
up between the system and environment then lead to 
an effectively irreversible, non-unitary dynamics of the 
reduced state of the quantum system. This is often clas- 
sified into two fundamental processes, the intuitive re- 
laxation of the system to thermal equilibrium, caused by 
energy exchange with the environment, and the destruc- 
tion of quantum mechanical coherence between classical 
system states, known as decoherence. 

In situations where the environment is weakly coupled 
to the system, methods such as Redfield or Lindblad the- 
ory can be applied^ and the dynamics of the reduced 
state can be described using simple time-local master 
equations. However, when the coupling is strong, or if 
the decay of environmental correlations is slow, then the 
Born-Markov approximation on which these techniques 
are based will fail, and a more sophisticated treatment of 
system-bath correlations and bath-memory effects is re- 
quired. The development of methods to treat such cases 
has recently been necessitated by interest in a wide range 
of quantum systems in which the environmental interac- 
tions and dynamics are non-trivial, leading to complex re- 
duced system state dynamics which are intermediate be- 
tween semi-classical energy relaxation and purely quan- 
tum coherent (wavelike) motion. Important examples 
include superconducting and spin quantum dot qubits 
for quantum computation^— and the recently discov- 



ered long-lasting wavelike motion of excitons in photo- 
synthetic pigment-protein complexes^— 

The extensively-studied spin-boson model has estab- 
lished itself as perhaps the most important system for 
developing theoretical concepts and numerical techniques 
through which we may understand the microscopic be- 
haviour of open quantum systems in all of the regimes 
mentioned above. The model consists of a quantum 
two-level system (TLS) that is coupled to an environ- 
ment made up of a bath (continuum) of harmonic oscil- 
lators. The environment and its couplings to the system 
are characterised by a spectral function J(uj), to be de- 
fined later. In many important applications, the spec- 
tral function behaves as a power-law at low frequencies 
J(uj) cx lu s , and this is often used to classify system- 
environment coupling types into three distinct groups: 
sub-Ohmic (s < 1), Ohmic (s — 1) and super-Ohmic 
(s > 1). Despite its apparent simplicity, the dynamics of 
the spin-boson problem cannot be solved exactly, and the 
extremely interesting, and as yet unexhausted, physics in 
the model continues to drive research into its properties 
and the potential implications for quantum devices and 
bio-complexes 

The super-Ohmic case is arguably the simplest to de- 
scribe and applies to a wide variety of physical sys- 
tems, such as in the electron-phonon interactions of 
impurities in solids-^ and quantum dots^ 3 . - — and in 
atom-photon interactions. 11 However, even in this case 
a crossover from coherent to incoherent dynamical be- 
haviour is expected as the environmental influence be- 
comes strong ) 22 i 27 ' 28 and the simplest weak system-bath 
coupling treatments will then fail. The sub-Ohmic 
case leads to strongly non-Markovian dynamics and 
also contains a quantum phase transition in the ground 
state^ir— Sophisticated numerical methods, capable of 
treating the many-body correlations which drive these 
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phenomena ; 32 ' 33 ' 39 — are often required to look at this 
case, and several artificial systems have been proposed 
in which these effects could be observed . 45 ' 46 The Ohmic 
case lies on the boundary of the non-Markovian, many- 
body physics of the sub-Ohmic regime and the effectively 
Markovian (though not necessarily weak-couplin g 22 ' 27 ' 28 ) 
physics which emerges in the super-Ohmic case. Its im- 
portance lies in the combination of its non-trivial dy- 
namical properties, which encapsulates a number of dif- 
ferent phases, and the many physical realisations of this 
type of environment . 1 ' 2 ' 1 1 A particularly topical example 
is found in exciton transport in pigment-protein com- 
plexes, where the (Ohmic) over-damped Brownian os- 
cillator model is widely employed as a simplified way to 
model solvent environments and protein fluctuations ! 16 ' 47 

To treat the dynamics of the spin-boson model, var- 
ious advanced numerical and analytical methods have 
been proposed and applied*^ Amongst these, a popu- 
lar and powerful approach has been the polaron trans- 
formation technique^ ' 22 ' 48 ' 49 which uses a unitary trans- 
formation of the spin-boson Hamiltonian to treat part of 
the system-bath coupling non-perturbatively, and then 
employs perturbation theory in the residual system-bath 
couplings to derive a dissipative master equation for the 
TLSi^ 7 .^ 8 -^— The perturbative treatment of residual in- 
teractions essentially drives thermalisation amongst the 
renormalised states of the non-perturbative part of the 
transformed Hamiltonian, and at low temperatures it is 
important that the transformation correctly captures the 
ground state of the system embedded in the environ- 
ment. Unfortunately, the naive application of standard 
polaron theory to Ohmic environments fails, as all co- 
herent matrix elements are renormalised to zero in the 
non-perturbative part of the Hamiltonian for all coupling 
strengths. This leads to purely incoherent dynamics in 
a time-local master equation approach^ even though 
more sophisticated numerical and analytical techniques, 
as well as common sense, show that damped coherent 
dynamics take place at sufficiently weak coupling^ 7 - An 
improvement on the standard polaron theory is the vari- 
ational polaron treatment of Silbey and Harris (SH) ^2r£2. 
in which the zeroth-order Hamiltonian may describe co- 
herent dynamics for weak coupling. However, for bi- 
ased TLSs this theory predicts an unphysical, discontinu- 
ous crossover to incoherent dynamics at a finite coupling 
strength^ 

The failure of these approaches lies in the improper 
choice of zeroth-order Hamiltonian, which at low tem- 
peratures results in the system relaxing to a ground 
state which is qualitatively different from the true ground 
state. In this paper we propose a new ground state 
ansatz for the biased Ohmic spin-boson model, and show 
that it predicts results in excellent agreement with non- 
perturbative treatments based on the Bethe-ansatz for 
this problem; 6 - 1 - As well as being conceptually simpler 
than these often costly techniques, our ansatz is itself 
based on a generalised polaron-type transformation that 
permits the nonequilibrium dynamics of the TLS to be 



explored using many of the advances recently made in 
polaron theory^^ 5 -^^ 4 - 

The paper is set out as follows. In Section [H] we in- 
troduce the spin-boson model and formulate the prob- 
lem. We then present our microscopic ansatz for the 
ground state and compare it to the SH theory and the 
exact Bethe-ansatz solutions of an equivalent theory. In 
Section |m] we demonstrate how the ansatz can be re- 
cast as a unitary transformation of the original problem, 
and derive the effective Hamiltonian with which we then 
compute the TLS dynamics. These results and the com- 
parison with the other theories are then discussed. In 
Section IIVI we briefly comment on the application of our 
ansatz to non-Ohmic spectral densities, before summaris- 
ing in Section [V] 

II. MODEL 

The spin-boson Hamiltonian can be written as 
H = |cr 2 - — a x + ^ ^k&k&k +P>y"] 3k (& k + fo k), (1) 

k k 

describing a TLS characterised by a bias e and tunneling 
amplitude A between basis states |1) and |0), linearly 
coupled to a harmonic oscillator bath of mode frequen- 
cies Wk, with strengths (assumed real). The standard 
Pauli operators used above are a z = |0)(0| — |1)(1| and 
<7 X = |1)(0| + |0)(1|. As has been well established by pre- 
vious studies) 1 -! 2 - for Gaussian initial states the effects of 
the environment on the reduced state of the TLS in this 
model are completely determined by the spectral func- 
tion J(oj), defined as J {to) — Sk^k^^k — w )- I n this 
paper we will only consider the Ohmic spectral density, 
which we parameterise as 

J{u)) = {a/2)u)6{u) c -u)), (2) 

with a being a dimensionless measure of the system- 
environment coupling strength, while 6(x) is the Heavi- 
side step function that provides a cut-off to the spectral 
function at a typical frequency of ui c . 

A. Ground state ansatz 

For computing non-perturbative dynamics in a stan- 
dard polaron-type theory, it is essential that the zeroth- 
order Hamiltonian possesses a ground state which is a 
very good approximation to the true ground state. We 
shall first construct such a ground state ansatz, and 
present an equivalent generalised polaron transforma- 
tion in Section IIIII The ansatz is composed from the 
basis {|0) JTk £, ( a k,o)|vac), |1) n k D ( a k,i)|vac)}, where 
| vac) denotes the vacuum state of the bosonic bath and 
D(ttk) = exp[ak(fo k — &k)] are bosonic displacement op- 
erators^ 5 - 
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In the new (restricted) basis, we may write the spin- 
boson Hamiltonian as 



H 



(e + A 



A 



* r ~ 



(3) 



where a z = |0) JIk D(a kj0 )|vac) (vac| ]J k D(-a k ,o)(0| - 
1 )rik D ( Q! k,i)|vac)(vac|nk D (- a k,i)(l|, 
Dk a k : o(^kak,o + 25 k ), Ai = J2k a k,i(wkak,i - 2g k ), 
and the tunneling term has been renormalised such that 

A r = A(vac| J^D(±(a k ,o - ak,i))l vac ) 



Aexp 



"oX^" 1 *' - «k,l) 



(4) 



At this point, the effective Hamiltonian of Eq. Q can 
be diagonalised and the ground state energy minimised 
as a function of ak,o and «k,i to find an optimal set 
of displacements in the ground state. For sub-Ohmic 
baths with s < 0.5, this procedure has been shown to 
correctly capture the physics of the mean-field quantum 
phase transition of this model, and reproduces results ob- 
tained by non-perturbative numerical techniques* 35 How- 
ever, for Ohmic systems we find that this procedure is 
unstable and deviates from the well-established results 
for the Ohmic case for a > 1/2. The cause for this is 
currently under investigation. Curiously, the correct be- 
haviour is found for the unbiased case (e = 0) by the vari- 
ational transformation of Silbey and Harris,— in which 
the constraint ak,o = — ctk.i is imposed. However, as will 
be discussed further in Section III Bt the SH procedure 
fails for the biased case. 

We now propose a non-variational ansatz for the mode 
displacements for the biased Ohmic spin-boson model, 
appropriate for the scaling-limit (A/uj c -C 1), which in- 
corporates features of both the sub-Ohmic and SH ground 
states, yet fixes the pathologies associated with both the- 
ories in this case. The proposed displacements are given 
by,^ 



"k,0 



ak,i 



.9k(e - ^k) 
Wk(wk + x) ' 
,9k(e + ^k) 
w k (w k + x)' 



(5) 
(0) 



where \ — vA 2 +e 2 . This leads to a self-consistent 
equation for the renormalised tunneling given by 



J(oj) 



A r = Acxp -2 / duj- , 

/o {u + xrJ 



which for J(ui) of Eq. becomes 



A r = A | — - — | e 

,X + ^c 



(7) 



(8) 



Eq. (J7J self-consistently predicts the renormalised tun- 
neling matrix element A r . It should be noted that in the 



presence of a bias, the integral is essentially cut off at low 
frequencies by the dynamical energy scale \, which is al- 
ways non-zero for e ^ 0. This means that in the biased 
case a finite solution to Eq. © can always be found, and 
A r is thus a continuous function of a. Similarly, we find 
(again for J(uj) of Eq. (0) 



Aq = I duj 
Jo 

auj c 



J(uj)(e-uj)(2x + e + uj) 
-(e 2 + X (2e- Wc )), 



J(<j)(e + uj)(e-u]-2x) 



Ai = I du 
Jo 

OtUJc 



2X(X + ^ C ) 



uj{u) + x) 2 
(e 2 - X (2e + w c )), 



(9) 



(10) 



which gives R = A — A\ = 2auj c e/ (% + w c ). Hence, we 
may now write the ground state energy in the Ohmic case 
as 

1 / auj c (e 2 - X uj c ) 



An = - 



x(x + ^ c ) 



n 



(ii) 



where rj = \/ A% + e 2 (l + 2auj c / (x + ^c)) 2 , while the 
ground state magnetisation, M — (cr z ), and coherence, 
(<x x ), become 

e(l + 2auj c /(x + u c )) 



M 



(12) 



and 



l*o) 



(or.) = VT^IP 

elf is i 

,o)|vac) 



(13) 



The approximate ground state itself is written simply as 
R + e-T] 



y/Al + (R + e- V y 
A r 



|l)n^K,i)|vac). 



(14) 



B. Comparison with the Silbey-Harris approach 

We shall now compare our ground state ansatz 
[Eq. ([H]) ] to that given by the variational treatment 
of Silbey and Harris.^ ~ 60 i 67 In the present formalism, 
the SH variational ground state is obtained by setting 
"k,o = -«k,i = -gk/(w k +A| H /xsH) inEqs. © and (gj), 
where Ash is the renormalise d tunnelin g element found 
in the SH theory, and %sh = i/Ag H + e 2 . These displace- 
ments can be obtained by minimising the ground state 
energy of the Hamiltonian of Eq. ([3]) with respect to ak,o 
and «k.i! subject to the constraint ak,o = — a k.i (which 
implies that R — )• 0). In doing so the self-consistent equa- 
tion for Ash becomes, 



^SH 



A exp 



duj 



J(w) 



(w + A 2 H / XSH ) 5 



(15) 
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FIG. 1. (a) Magnetisation, M, as a function of system-bath 
coupling, a, for the Ohmic spin-boson ground state. Re- 
sults from the Silbey-Harris variational treatment are shown 
as open circles, while the solid curves are plotted using the 
ansatz presented in this work, see Eqs. (|12fl and (|14[l . Part 
(b) shows the difference in the ground state energies predicted 
by the Silbey-Harris theory and the present ansatz (both of 
which are negative). In both plots A = 10~ 2 a; c , with the 
arrows indicating values of e = 0.5A (blue), e = 0.1A (red), 
e = 10" 2 A (green), e = 10" 4 A (orange), and e = 10" 8 A 
(purple) in decreasing order. 



which leads to 
A SH = A 2 5 exp -2— , (16) 

V A§ H + XSH^c J L Ag H + XSH^c J 

again in the Ohmic case. Note that the low energy 
cut-off scale is now given by Ag H /x SH which can self- 
consistently vanish above a critical coupling strength. 
This is the essence of the SH theory at strong coupling. 
The SH ground state energy is given by 

, If otXsnUc \ ri ~s 

2 V(Ash +XSHW C ) ) 

while we also find Msh = (c z )sh = — c/xsh and 
(Osh = A| h /xshA. 

In Fig. Q] (a) we plot the magnetisation of the present 
spin-boson ground state ansatz (from Eq. (TTiZj)) as a func- 
tion of the dimensionless coupling a, and compare with 
that given by the SH method, for various values of the 
bias e. As is now well known, for small e the SH magneti- 
sation, Msh, displays an unphysical discontinuous 'jump' 



to Msh = — 1 (corresponding to a fully localised ground 
state) for some value of a < 1. This behaviour can in 
turn be attributed to a discontinuous change from non- 
zero Ash to Ash = in the SH theory as a is var- 
ied. From Fig. [T] (a), however, we see that the ansatz 
presented in this work leads to no such 'jumps' in the 
behaviour of M, and the magnetisation smoothly ap- 
proaches — 1 with increasing a. This behaviour is in 
agreement with that found from various advanced numer- 
ical methodS ) 32 ' 33 ' 68 i 69 as we shall show explicitly below 
using results derived from the Bethe-ansatz. 

In Fig. [1] (b) we plot the difference in ground state 
energies predicted by the two theories, Ash — Ao, as a 
function of a, for the same values of the bias as in part 
(a). We see that for all parameters considered here, the 
present ansatz corresponds to a ground state approxima- 
tion with lower energy than that given by the SH theory 
(both Ash and Ao are negative for these parameters), 
suggesting that the state given by our ansatz is indeed 
a better approximation to the true ground state. Inter- 
estingly, we see that the difference in the two ground 
state energies is maximised through the crossover region, 
where Msh is changing most rapidly. We also note that 
the unusual behaviour of Ash — Ao near the peak differ- 
ence is due to the discontinuities present in the values of 
Msh and Ash as a is increased. 



C. Comparison with the Bethe-ansatz 

As discussed in Refs. HI. 13. loll. lo8l l7fjl - [73l a mapping ex- 
ists in the scaling limit (A/uj c <C 1) between the Ohmic 
spin-boson model, the anisotropic Kondo model, and a 
range of interacting resonance-level models. Exploiting 
this mapping and the existence of exact Bethe-ansatz so- 
lutions for the resonance-level model, explicit formulae 
for the properties of the biased spin-boson model can be 
obtained. To compare to the present ansatz results, we 
first take the Bethe-ansatz expressions for (o~ z ) given in 
Eq. (CI) of Ref.H 

In Fig. [2] we plot a comparison of the magnetisation 
predicted by our ansatz (solid curves) , and that given by 
the Bethe-ansatz (crosses), for the same parameters as 
in Fig. Q] We see that the agreement between the two 
methods is generally very good, though there are small 
deviations in the crossover behaviour at the weakest bi- 
ases in plot (a). However, if we move further into the 
scaling- limit, by reducing the ratio A/w c as in plot (b), 
we find excellent agreement across all bias values. This 
improvement with reducing ratio of A/w c is not too sur- 
prising, as we expect both the present ground state and 
the Bethe-ansatz to be most valid in the scaling-limit. In 
the former case, due to the restricted basis from which the 
present ground state is constructed, in the latter case due 
to the mapping that is utilised from the spin-boson model 
to the anisotropic Kondo model, which is only valid in 
the limit that A/w c becomes very smalh 61 ' 70 " — Still, it 
is remarkable how well the two methods agree, particu- 
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FIG. 2. Magnetisation, M, as a function of system-bath cou- 
pling, a, calculated from the present ansatz (solid curves) and 
calculated using the exact scaling-limit expressions obtained 
via the Bethe ansatz (crosses), presented in Ref. loll The ar- 
rows indicate decreasing values of e/A, as in Fig. [1] The two 
plots correspond to different values of A/cu c , as indicated. 
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FIG. 3. Ground state spin coherence (a x ) as a function of 
system-bath coupling, a, calculated from the present ansatz 
(solid curves) and calculated using the exact scaling-limit ex- 
pressions obtained via the Bethe ansatz (crosses), presented 
in Ref. [611. The arrows indicate the bias values e = 0.5A 
(blue), e = 0.1A (red), and e = 10" 2 A ( green), in decreasing 
order. The two plots correspond to different values of A/w c , 
as indicated. 



larly over the sharp crossover region, given the simplicity 
of our proposed ground state [Eq. (|14|)] in comparison to 
the full Bethe-ansatz. 

Similar behaviour is seen for the ground state spin co- 
herence {(J x ), as shown in Fig. [31 Here, our ansatz leads 
to good qualitative agreement with the Bethe-ansatz but 
consistently predicts slightly lower coherences. Again, 
these deviations become smaller as we move further into 
the scaling limit. We note that the slightly weaker agree- 
ment for these quantities is likely to be due to the fact 
that they are non-universal and vanish in the scaling 
limit.— The finite values we obtain thus depend on details 
of the high frequency cut-off procedure. While the Bethe- 
ansatz results for the Kondo/resonance level model are 
still exact, the correspondence between them and the 
spin-boson model results is also dependent on the details 
of the high-frequency regularisation used in the mapping 
that links them in the scaling limit. As these details do 
not necessarily coincide for our ansatz and the Kondo 
mapping, numerical differences in coherences are likely 
to occur. However, as shown below, the leading-order 
functional forms of the non-universal properties are the 
same, and the numerical agreement for the magnetization 
(a universal property) is extremely good. 

We can also compare the ground state energy predicted 



by our ansatz [Eq. (fTTj) ] with that given by the Bethe- 
ansatz (see Eqs. (C3)-(C9) of Ref. Wu \. Examples are 
shown in Fig. [H where we need to add a term equal to 
aoj c /2 to Eq. (fTTj) to be consistent with the ground state 
energy definition in Refs . Ill [fill 1701473 . Once more, we see 
that the agreement improves as the ratio A/w c decreases. 

1. Analytical results 

Figs. [2j 02 and [4] suggest that, in the scaling-limit at 
least, the present ansatz provides a very good approxi- 
mation to the true ground state of the Ohmic spin-boson 
Hamiltonian. We now investigate this further by com- 
paring the analytic expressions obtained from this work 
with those presented in Ref. |6lj. To recap, the magneti- 
sation of the ground state in the present theory takes the 
form 

e(I + 2aw c /( X + w c )) 



where r\ = ^/A 2 + e 2 (I + 2auj c / (x + ^c)) 2 and \ — 
A 2 + g 2 . In the scaling limit, x/w c — > 0, simple forms 
for the magnetisation and other parameters such as A r 
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FIG. 4. Ground state energy, A, as a function of system- 
bath coupling, a, calculated from the present ansatz (solid 
curves) and calculated using expressions obtained via the 
Bethe-ansatz (crosses), presented in Ref. loll The arrows in- 
dicate the bias values e = 0.5 A (blue), e = 0.1 A (red), and 
e = 10 _2 A (green), in decreasing order. The two plots corre- 
spond to different values of A/cj c , as indicated. 



can be obtained. When e <C A r , we obtain to lowest 
order in e/u; c , 



M 



e(l + 2a) 



(18) 



The self-consistent equation for A r in the scaling limit 
for weak bias is 



Ae« | ^ 

UJ C 



which can be solved analytically and substituted into Eq. 
(UHl to get, 



M oc 



e 

Tk 



(19) 



where we have used the definition of the Kondo tem- 
perature, T K = A (A/D) a ^ 1 ^ a , and the further relation 
between D and w c given in Eq. (C8) of Ref. [6l| to facili- 
tate an easier comparison with the Bethe-ansatz results. 
In this limit, we also obtain for the coherences, 



/ \ ~ Tk / i i a 



(20) 



For strong bias (e ^ Tk), we find that 
A. r w A 

and the magnetisation and coherence become 
M+l cc (— V (— 2 '~ 2 



and 



2a-l 



(21) 



(22) 



respectively, where, by our choice of notation, M = — 1 
is the magnetisation of the TLS in the absence of bath 
coupling in the limit e/A»l. In both limits, these func- 
tional forms for the observables of the TLS coincide with 
the predictions of the Bethe-ansatz, and the cumbersome 
numerical prefactors we have omitted appear from our 
comparisons in Figs. [HO and 2] to be very close as well. 



III. UNITARY TRANSFORMATION 
APPROACH 

Having now shown that l^o) given in Eq. (fl"4")l can 
provide an excellent approximation to the ground state 
properties of the Ohmic spin-boson model, let's return 
to the original spin-boson Hamiltonian of Eq. ([T]) to ex- 
plore a unitary transformation approach to the problem. 
We shall show how this can be made equivalent to the 
ground state ansatz method outlined above, and how the 
transformation also provides a basis for computing the 
dynamics of the TLS beyond weak system-bath coupling. 

We perform the unitary transformation, H' = 
e s He~ s . where 



jjD(±* kj0 )+ii>(iin £, ( ± ^.i)' ( 23 ) 



to give 



H' = *--(T z - ^-cr x + ^Ldkblbk - ^-(a x B x + cr y B y ) 



+ \0)(0\(B zfi + A' o )-\l)(l\(B zA -A' 1 ). (24) 



Here, 



B , 



B+ + B_- 2B 



, By 



B_ - S_| 
2i 



and 



B z ,o = y^(ffk ~ Wk<5k,o)(&ic + K), 



B z s = J^(.9k + ^A,i)(6 k + ht), 



7 



while 



A' = ^ <5k,o(wk<5k,o - 2gk), 



k 



<$k,l( w k<5k,l + 25k), 



and the renormalised tunneling is now given by 
A r = BA, where we have defined B = trs(_B±ps), 



terms of B± 



n k ^(±(^o 



4,0) 



and a thermal equilibrium bath state ps — 
exp[-/3 Ek ^k&I&k] I tr B (exp[-/3 £ k w k 6 k 6 k ] ) . The 
inverse temperature is /3 = l/fe B T. 

Immediately, we see that if we choose <5k,o = — Cfck,o 
and (5k. 1 = -«k.i, we get the same definitions of Aq and 
Ai as those used previously (i.e. A' Q — > A , A[ — > A\). 
Hence, we take 



•>k,0 



?k,l 



-«k,0 



-«k,l 



gk(^k - e) 
Wk(w k + X) ' 
g k (^k + e) 
w k (w k + x) ' 



(25) 
(26) 



where x = -y/A 2 + e 2 as before. With these definitions, 
the renormalised tunneling becomes 



A r = A exp 



A exp 



£>k,i) 2 coth 




: , J w , fa 

dhj- coth — 

(w + x) 2 2 



(27) 



in the continuum limit, which is a finite temperature gen- 
eralisation of Eq. (O . 

Let us now split the transformed Hamiltonian as W = 
H'q + H'j, where 



and 



A, 



k 



w k & k &k+|0)(0|A) + |l)(l|Ai, 

(28) 



A 



^1 — — — (fx-Ba; + CTyBy) " 



|0)<0|B,, -|l)<l|B,,i, (29) 



which ensures that (ifj)^ = tr[i?je~' 3ff o/ tr ( e ~' 3ff ")] = 
0. The Feynman-Bogolioubov upper bound on the 
free energy^ A B = -(1//3) ln[tr( e -^ H o)] + {H'j)^ + 
O((Hj 2 ) H ' ), then becomes 



A + A 1 _ 1 
2 



In 



2 cosh 



2 



(30) 



where 77 = ^A 2 + (i? + e) 2 , with i? = A n — Ai, exactly 
as before, and we ignore contributions from the free bath 
Hamiltonian as we are interested only in the free en- 
ergy of the TLS. In the zero temperature limit (/? — > 00) 
Eq. (J3DJ) becomes 



A, 



An 



-(A + A 1 -rf), 



(31) 



in agreement with the ground state energy of Eq. (ITT1) . 

Furthermore, approximating the thermal state density 
operator in the transformed frame as ptn = e~@ H °/Z, 
where Z — ti(e~^ H °), we find that the thermal expecta- 
tion values of the Pauli spin operators may be written 



(di) = tr (e s cr i e s p th ) 
Hence, we find at finite temperature 

(^) = -^±i) t anh^ ; , 



and 



(ctj) = — 7- tanh ( 



(32) 



(33) 



(34) 



both of which agree with our previous expressions 
(Eqs. (fT2"j) and ([IB"]) , respectively) in the zero tempera- 
ture limit. 



A. Dynamics - Master equation derivation 

The advantage of the unitary transformation approach 
is that we are now in a position to efficiently explore the 
reduced TLS dynamics within the model, by deriving a 
master equation in the transformed representation. The 
general philosophy is that given we have shown that our 
zcroth-order Hamiltonian H' a provides a good approx- 
imation to the model ground state over a wide range 
of system-bath coupling strengths, we should expect (at 
low temperatures at least) that the effects of the pertur- 
bation H'j remain small over this parameter range too. 
Thus, we shall treat Hj in low-order perturbation the- 
ory within the scaling limit, and derive a second-order 
master equation that should be valid at strong system- 
environment coupling as well as in the more usual weak- 
coupling regime. 

Using our partitioning of the total Hamiltonian 
given by Eqs. (|2"5)l and (f2T)|) . together with the time 
convolution-less projection operator technique^ we ob- 
tain an interaction picture time-local master equation of 
the form 



dp v (t) 
dt 



drtTBdH'j^^H'jit-^^^pB}}) 



(35) 



valid to second order in Hj. Here, p v (t) is the re- 
duced density operator of the TLS degrees of freedom 
in the transformed frame interaction picture, and the 



iH' n t 



interaction Hamiltonian is Hj(t) = e ot H'je 
- (A/2) (oi (t)B x (t) + o y (t)B y (t)) + a z (t)B z (t) + tBj (t) , 
where we have defined the bath operators B z — 
(l/2)(B, i0 +fl,,i)andBj - {\/2){B zfi -B zA ). We write 
ai (t) = e iH sVie- iH s* and B,(t) = e iHBt B t e- iHBt , for 
i = I,x,y,z (where 07 = 1), with H' s = ((e + R)/2)a z — 
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(A r /2)(T X and Hb = Ek w k^k- In deriving Eq. (1351) . 
both the reference state in the projection operator and 
the environment initial state (in the transformed frame) 
are chosen to be the thermal equilibrium state p B , result- 
ing in the absence of any inhomogeneous termsi 11 ' 52 ' 53 ! 63 
We now insert the interaction Hamiltonian into 
Eq. (|35[) . take the trace over the environment, and make a 
secular approximation to remove fast oscillating terms ^ 
This last simplification is made in order to ensure that 
our system relaxes to its ground state in the long time 
limit (at zero temperature), which is appropriate given 
the close agreement we find between the properties of 
our ansatz ground state and those of the Bethe ansatz. 
Moving back to the Schrodinger picture, we then find 

Pv = - 2 I( e + R ) a z - A r a x ,p v ] 

-^E{^4,A^] r -(V) 

A 



[^,A,4 iA p,]r yy (A,<)+H.c.} 



A 



2 E { A l\PA - [A z ^Al xPv ])T yz {\, t) 

A 

+ S x .o[A y ,x,p v }T yI (X,t) +R.c.} 

^2{[A x , x ,Al tXPv ]r„(x,t) 



5x.o[A z ^,p v ]r zI (X,t) + U.c.}, 



(36) 



where Sx.o is the Kroneckcr delta, and we have decom- 
posed the system operators as Oi(i) = Yl\ e%Xt Ai,\, for 
A = 0, ±rj. In terms of 6 = arctan(^^), we have A x $ — 



sin (sin 9a x + cos0a z ), A x<n 



cos 
2 



-(cos 9a x 



sin 9<j z ), A yfi 



0, A. 



| (cos 9a x 



io y — sin 8<j z ), 



A Z fi — cos9($m.9a x +cos6a z ), A z , v = s §^(— cos#ct, 
i(j y + sin 6<7 z ), while — 1 

functions 



AT The bath response 



dre iXT 



Ay(r), 



(37) 



are defined in terms of the correlation functions Ay (r) = 
tr B(Bi(r)Bj(0)p B ), for i,j =I,x,y,z. 

We proceed by writing the transformation parameters 

as 

5k 



?k,o 



3k. i 



Ck,0 : 

, 9k 

Sk,l 5 



(38) 
(39) 



such that Ck,o = (wk - e)/(w k + x) and Ck,i = -(wk + 
e)/(wk + x) f° r our ansatz [see Eqs. (f25j) and ([2"6)) ]. We 
may now re-write the renormalised tunneling as 



A, 



A exp 



'^■ilCk.o-Ck^coth^ 
k -k 1 , 



2 4^ w 



1 



(40) 



where C±M = CoM±Ci(u;), with CoM = (w-e)/(w+x) 
and Ci(w) = ~ ( w + e ) /( w + x) i n the continuum limit. 
In terms of C±( w ) the correlation functions are found to 
read 



A xx (t) 
A to( t ) 



(C(t)+G(t)-2B 2 ), (41) 
(C(r)-G(r)), (42) 



where 
C(r) = exp 

G(t) — exp 



x ((1 — cos wr) coth fiw/2 + i sin wt) 



x ((1 + cos wr) coth/3a;/2 — i sinwr) 



(43) 



and -B = A r /A, while 



A, z (r) = -B /°° dw^C-H fl 
Jo w V 



C-M 



/>oo 

A 2//( T ) =+ B duJ 
JO 



x (sin lot coth /3w /2 + i cos wr) , 

JH C-HC+H 



x (sin ojt coth /3w /2 + i cos wr) , 

C-H 



A zz (r) = / d W J(u;) 1 



x (cos wt coth (3ui /2 — i sin wr) , 



x (cos wt coth /3w /2 — i sin wr) , 

(44) 

and we have already used the fact that A xy {r) — 
A yx (r) = A xz (t) = A zx (t) = A xJ (r) = A Ix {t) = 
in deriving Eq. ([36]). We note that Eq. (|36| and the 
subsequent definitions are quite general, and the (secu- 
lar) master equation corresponding to any polaron-likc 
transformation on the Hamiltonian can be obtained by 
the appropriate choice of the functions (o(uj) an d Ci( w )- 
For example, the master equation corresponding to the 
SH transformatio n 62 ! 63 is given by setting (q H (u) = 
u/(u} + A sh /xsh) and Cf^M = + A sh /xsh), 

which we shall use for comparison to our ansatz below. 

We now solve the ansatz master equation numerically 
to determine the dissipative spin dynamics, concentrat- 
ing on the zero temperature limit. At zero temperature, 
and for an Ohmic spectral density, the correlation func- 
tions defined above do not decay on a rapid timescale 
compared to the system dynamics we wish to capture. 
As such, in the following we cannot extend the integra- 
tion limits in Eq. (|37[) to infinity. We consider the system 
to be initialised in state |0) at time t = 0. 
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In Fig. [5] we plot the population dynamics, (<r z ) t — 
tr(a z p v (t)), using (a) our present ansatz, and (b) the SH 
theory, as a function of the relevant scaled time, i.e. A r t 
and Ash^, respectively. The solid curves in (a) show 
coupling strengths ranging from a = 0.1 to a = 0.6 for 
the present ansatz, while in (b) the range a = 0.1 to 
a = 0.5 is plotted using the SH theory (in this case the 
a = 0.6 curve is not plotted as Ash = already for 
the SH theory). The other parameters are A = 10 _2 u; c , 
e = 10~ 2 A, and T = 0, which correspond to the central 
green curve in Fig. [T] For the smallest coupling strength 
(a = 0.1), we see that the dynamics calculated from the 
present ansatz and from the SH transformation agree al- 
most perfectly. This is unsurprising, as we know from 
our previous comparisons that for these parameters the 
ground states given by the two methods are very similar. 
Nevertheless, discrepancies do begin to become apparent 
between the two methods at larger a, as shown by the 
curves corresponding to a — 0.25, a = 0.4, and a = 0.5. 
In particular, we see that in these cases the ansatz steady- 
states settle at long times to slightly lower values of (a z ) 
than the SH dynamics predicts. This is to be expected, 
given the corresponding ground state magnetisation be- 
haviour shown in Fig. Q] and the fact that the two-state 
system relaxes towards the relevant ground state at long 
times for T = in the secular master equation. 

At coupling strengths beyond a ~ 0.55, the SH the- 
ory predicts a complete renormalisation of the tunnel- 
ing strength to zero, i.e. Ash — > 0, while for a < 1 
the tunneling always remains non-zero when using the 
present ansatz. Hence, significant differences emerge in 
the population dynamics, with the two approaches dis- 
agreeing on both the population decay rates and steady- 
states. In particular, as a consequence of Ash — ► for 
a > 0.55, the SH theory always incorrectly predicts a 
system steady-state that is completely localised in this 
regime, which is not the case for our ansatz. For bi- 
ased systems at T = 0, this failing is also well-known 
to be true of the non-interacting blip approximation (for 
all a < l). 1,2 Again, the problem can be traced back to 
the incorrect choice of zeroth-order Hamiltonian in the 
SH theory (and in effect within the non-interacting blip 
approximation as well). 



1. The Toulouse point 

We have shown that our unitary transformation ap- 
proach (correctly) leads to significant changes in the spin 
dynamics at strong coupling compared to the SH theory. 
However, to assess the accuracy of the new dynamics, we 
now compare these results to the exact expressions that 
are known at the Toulouse point of the spin-boson model, 
which occurs at a — 0.5. At this coupling, the spin-boson 
model becomes equivalent to a non-interacting resonance 
level model in the scaling limit, and the spin dynamics 
can be solved exactly.^ 2 - 

The dotted curves in Fig. [5] show the exact analytical 



a) 1.0 
0.5 

<cr z >, 0.0 
-0.5 
-1.0 



1 

















5 10 15 20 25 

b) 1.0, v _ 

0.5 \ 

«r z > t 0.0 \mN — 1 / ~\ ~- ,< — ; — 

-o.5 vfeuiZ _jl ; 



5 10 15 20 25 

A S Ht 

FIG. 5. Magentisation dynamics (o~z)t computed within the 
secular approximation as a function of (a) A r t for the present 
theory and (b) Asnt for Silbey-Harris theory. Solid curves: 
the arrows indicate the coupling strengths a = 0.1 (blue), 
a = 0.25 (red), a = 0.4 (green), a — 0.5 (orange), and a = 0.6 
(purple, shown for our ansatz only) in increasing order. The 
dotted curves correspond to the exact (scaling limit) solution 
at a — 0.5, extracted from a mapping to the Toulouse Hamil- 
tonian, and computed from the formula given in Ref.0. The 
dashed curves show the same exact solution, but now plot- 
ted in units of Tkt, where Tk = A 2 /iv c is the Kondo scale 
for the Toulouse problem.— Other parameters: A = 10 _2 u; c , 
e = 10~ 2 A, and T = 0. 



solution for a = 0.5 taken from Ref. [1, as a funnction 
of scaled time A r t in (a) and Ash^ in (b). We see that 
the final value of (<j z ) within our ansatz coincides with 
that of the exact Toulouse point solution in (a), while, 
as expected, the SH theory dynamics in (b) tends to a 
substantially less magnetised, and incorrect, final state. 
We observe, however, that the dynamics of the Toulouse 
point solution appears to be slower than that predicted 
by our master equation approach. We suggest that the 
origin of these different decay rates could be related to 
the non-universal factors which appear in the ansatz pre- 
dictions for A r . Specifically, the exponential factor that 
appears in the self-consistent equation for A r [Eq. (jHJ)] is 
not a universal feature, in the sense that it depends on the 
way that integrals over the bath spectral function are cut- 
off for uj > uj c . For instance, using an exponential cut-off 
in the Ohmic spectral function, J{lo) — (a/2)uje~ UJ / u ' c , 
yields the following self-consistent equation for A r in the 
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scaling limit: 

A r = A e Q ( 1+ ^, (45) 

where 7 0.577 is the Euler-Mascheroni constant, 
as compared to A r = A(x/w c ) Q e Q for J(w) = 
(a/2)u>8(u> c — lo). In the absence of a bias, we see 
that the solution of the self-consistent equations then 
leads to an effective Kondo temperature Tk = A r = 
A [Af /ui c ) a ^ 1 ~ a \ where / is numerical factor (« e) 
which depends on the form of the cut-off. The factor / is 
non-universal, and this is also true of SH theory. In the 
exact solution presented in Ref. [||for the Toulouse point, 
the Kondo temperature that appears has / = 1, and is 
therefore smaller than the predicted by our ansatz or 
SH theory. Therefore, when we plot dynamics in units of 
A r t or Ash^, computed from Eqs. ([5]) and (IT5)l . respec- 
tively, we are not using the natural units for the Toulouse 
point solution. 

The dashed curves in Fig. [5] show the Toulouse point 
solution when plotted as a function of A 2 t/uj c , which is 
the appropriate Kondo temperature at a — 0.5. The 
ansatz and Toulouse point curves in (a) now show much 
closer agreement, whereas SH theory still captures nei- 
ther the right dynamical timescale nor, of course, the 
final magnetisation. 

These results suggest more generally that variational or 
ansatz-based polaron methods of the kind explored here 
may suffer from artefacts arising from the treatment of 
low frequency modes, which lead to dependencies of the 
effective Hamiltonian parameters on the details of the 
high frequency cut-off procedure. As far as we are aware, 
this has not previously been discussed or explored in de- 
tail, and is usually neglected. We have seen here, how- 
ever, that this may have a material effect on the dynamics 
of the system, and would therefore be an interesting issue 
to explore in future work. 

IV. NON-OHMIC SPECTRAL DENSITIES 

In this section we briefly comment on the application 
of our generalised polaron theory to non-Ohmic environ- 
ments. However, for reasons to be set out below, we be- 
lieve that the particular form of our ansatz is likely only 
to be useful in the Ohmic case. Although the sudden col- 
lapse of coherent tunnelling is not an issue with super- 
Ohmic baths, we could expect that our ansatz/unitary 
transformation may still lead to different dynamics in 
such a case compared to the SH theory, for example 
through the change in dynamical cut-off in the renor- 
malised tunnelling matrix element and the bath-induced 
bias (which does not feature in SH theory at all). How- 
ever, repeating the calculations of Section III Al for spec- 
tral densities of the form J(w) cx oj s (with s > 1) shows 
that the effects of asymmetric displacements in the renor- 
malised tunnelling amplitude and bath-induced magneti- 
sation vanish in the scaling limit. As renormalisation 



effects converge to those predicted by adiabatic renor- 
malisation, or, equivalently, standard polaron theory, (as 
they also do in SH theory) we therefore expect this the- 
ory to coincide with full polaron theory for supcr-Ohmic 
baths as uj c — > 00. 

As we have previously mentioned, the asymmetric dis- 
placement terms are extremely important in the sub- 
Ohmic case, and may appear spontaneously even with- 
out a bias at strong coupling and T = Repeating 
the calculations of Section III Al shows this time that the 
effective bias induced by the environment is proportional 
to (cj c /x) 1_s , and thus diverges in the scaling limit. 
This would predict complete localisation for any alpha 
and would be inconsistent with the known existence of 
a quantum phase transition in sub-Ohmic systems^ - — 
However, the variational ansatz proposed by Chin et al. 
in Ref. 35 regulates this pathology, and could also be used 
as the basis for a variational unitary transformation that 
is similar to the dynamical approach pursued here. This 
will be discussed in a forthcoming work. 



V. CONCLUSIONS 

We have presented a new ansatz for the ground state 
of the biased Ohmic spin-boson model which is of lower 
energy than the Silbey-Harris state and cures the prob- 
lem of the discontinuous collapse of the coherent tunnel- 
ing matrix element as the coupling strength to the en- 
vironment is increased. The key differences between our 
ansatz and the SH variational state result from the dif- 
ferent forms of the displacements «k,o an d Qfk.i which are 
taken in these ground states. Firstly, the SH ground state 
has strictly (anti) symmetric displacements and the bias 
only appears in the low frequency energy scale Ag H /%sH- 
In our ansatz the displacements are asymmetric due to 
the inclusion of the bias in the numerators of Eqs. ([5]) and 
([6]) . Similar asymmetric displacements were shown to be 
essential for describing the spontaneous magnetization of 
the unbiased sub-Ohmic spin-boson model which charac- 
terises its quantum phase transition, and in effect causes 
an extra bath-induced bias to be seen by the TLSj 35 i 38 
Moreover, we have shown explicitly that asymmetric dis- 
placements lead to the appearance of the correct dynam- 
ical low energy cut-off (x) in the self-consistent equation 
for A r for the biased Ohmic case. The close agreement 
of results obtained through our ansatz and the Bethe- 
ansatz evidences that our ansatz is capturing a consid- 
erable amount of the essential physics of the problem. 
However, the detailed links between the rather simple 
form of the spin-boson model ground state we have given 
and that of the equivalent Kondo model remain to be an- 
alyzed, particularly for non-universal observables, such as 
the spin coherence, and outside of the scaling limit. 

Importantly, by showing that this ground state can 
be obtained via a unitary transformation technique, it is 
now possible to treat the dynamics of the Ohmic spin- 
boson model using master equation techniques from po- 
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laron theory. The examples given in this article al- 
ready illustrate the dramatic differences in relaxation 
behaviour, particularly in terms of steady states, when 
the zeroth-order Hamiltonian correctly preserves a coher- 
ent tunneling matrix elements at all coupling strengths. 
This may be of relevance for parameter regimes found in 
some biomolecular complexes, where biases, environmen- 
tal coupling strengths and (bare) tunneling amplitudes 
are comparable^ conditions under which the differences 
between SH theory and the present theory are greatest 
(at T = 0). Although some deviation was found be- 
tween our results and the nominally exact dynamics at 
the Toulouse point, it is clear that our theory still per- 
forms better than SH theory at this point, and we were 
able to understand where these differences arise. While 
our ansatz cures one of the most important problems 
associated with SH theory in biased Ohmic systems, in- 
sights such as those given above also point the way to 



future refinements of the present theory. 

Finally, our ansatz could also be applied using the more 
advanced, non-Markovian master equation techniques, 
such as those presented and analysed in Refs. HH, [HI, 
[H, and [2!, which can also include bath relaxation ef- 
fects. These improvements in the handling of the dynam- 
ics could go some way to reducing some of the differences 
we found with the Toulouse solution, and along with ef- 
fects due to finite temperatures and moving beyond the 
scaling limit, are interesting areas for investigation in the 
near future. 
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